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abstract 


Computer simulations of realistic ion channel structures have always been challenging 
and a subject of rigorous study. Simulations based on continuum electrostatics have 
proven to be computationally cheap and reasonably accurate in predicting a channel’s 
behavior. In this paper we discuss the use of a device simulator, SILVACO, to build a 
solid-state model for KcsA channel and study its steady-state response. SILVACO is a 
well-established program, typically used by electrical engineers to simulate the process 
flow and electrical character- istics of solid-state devices. By employing this simulation 
program, we have presented an alternative computing platform for performing ion 
channel simulations, besides the known methods of writing codes in programming 
languages. With the ease of varying the differ- ent parameters in the channel's 
vestibule and the ability of incorporating surface charges, we have shown the wide- 
ranging possibilities of using a device simulator for ion channel simulations. Our 
simulated results closely agree with the experimental data, validating our model. 
© 2006 Elsevier Ireland Ltd. All rights reserved. 


1, Introduction 


residues on the wall, the surrounding ions, and the charges 
induced in the membrane wall by these ions [3]. 


Biological ion channels are transmembrane proteins that 
reg- ulate the ion flux through a cell, and thereby play an 
inte- gral role in cellular signaling mechanisms [1]. 
Modeling of ion channels has a long history, but the lack of 
detailed structural knowledge hampered the progress in 
this field. With advance- ments in molecular biology and X- 
ray diffraction techniques, a lot is now known regarding the 
structural details of ion chan- nels [2]. 

lons in solution execute random Brownian motion, 
and if an electric field is applied, they acquire drift velocity. 
It is the interaction between the ions and the electric field 
in the channel that decides the salient features of an ion 
channel. The effective electric field is a complicated 
interplay from var- ious sources including the membrane 
potential, the charge 


* Corresponding authors. Tel.: +1 610 758 4421. 


Among the varied approaches available for ion channel 
modeling, the most accurate representation of the system 
is provided by molecular dynamics (MD) simulations, 
where all the atoms in a system are treated explicitly [4,5]. 
However, MD simulations are presently limited to 
simulation times of 
not more than a few nanoseconds, due to the small 
timesteps (~femtoseconds) required to resolve individual 
ion trajecto- ries. Any reliable estimates of steady-state 
channel currents 
require simulation time durations of the order of 
milliseconds, thus preventing an experimental validation of 
MD simulations [6]. 

An alternative method is to follow the motion of individ- 
ual ions using Brownian dynamics simulations [7,8]. Here, 
it is assumed that the protein structure is fixed and 
water 
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molecules are treated as a continuum. Collisions between 
the ions and the surrounding water molecules are 
mimicked by random fluctuating forces plus an average 
frictional force. A major advantage of Brownian dynamic 
approach is the ability to allow direct simulation of ion-ion 
interaction and calculate an ion’s interaction energy with 
respect to the system during each time step [7]. But this 
has always been an involved, time- consuming process [9]. 

Using a mean _ field approximation, continuum 
electrostat- ics provides a rather simplistic alternative to ion 
channel mod- eling [10]. It involves solving the Poisson- 
Nernst-Planck (PNP) equations for the charge distribution 
inside the channel, and has proven to be very 
computationally cheap. Even though the PNP method uses 
a major approximation of reducing the ion-ion interactions 
to an interaction between an ion and a mean field, its 
results have shown to be in good agreement with the 
experimental data [11,12]. 

In this paper we _ present a novel method of 
simulating ion channels based on continuum electrostatics. 
We have attempted to represent the structure of a KcsA 
channel as a solid-state device, with its intrinsic material 
properties similar to that in the realistic case [13,14]. 
Simulations are performed on a_ solid-state device 
simulator, SILVACO International, pro- viding us with the 
steady-state behavior of the KcsA channel. The simulation 
results agree well with the available experi- mental data. 


2. Methods 


2.1. Poisson-Nernst-Planck (PNP) theory 

The PNP approach to ion permeation in membrane 
channels has been studied in numerous works [10,11,15]. 
The PNP equa- tions are similar to the drift-diffusion 
equations that describe the carrier dynamics in 
semiconductors [16]. The flux / of each ion species is 
described by the Nernst-Planck (NP) equation, which 
combines the diffusion due to a concentration gradient 
with the potential gradient [16,17]: 


j= ee VE (1) 
kT 


where D is the diffusion coefficient, z the valency of the 
ion, e the electronic charge, k the Boltzmann constant, 7 
the tem- perature, n the concentration of the ion, and ° is 
the potential. Also, the Poisson’s equation is written as 
{16,17]: 


EoVlE(r)V¢(r)] = - zen - Dex 


(2) 


ions 


where é(r) is the dielectric constant at an axial distance, r 
the sum over all ions gives the total mobile charge, & the 
dielectric constant of free space, and pex represents all the 
external fixed and induced charges. It is to be noted that, 
realistically, all the quantities in Eqs. (1) and (2) are 
varying along the ion chan- nel axis [12,13]. This makes it 
notoriously difficult to obtain any analytical solution of the 
PNP equations, besides in some very specific cases [10]. 
However, it is possible to solve Eqs. 


channel axis [17]. Even though this paper deals with 1D 
PNP computations, it has been shown that the PNP 
equations can be extended to 3D with a more realistic 
channel representa- tion and even closer fits to the 
experimental data [12]. 
2.2. The SILVACO simulator 
SILVACO is a solid-state process and device simulator that 
has the general capabilities for numerical, physically 
based, 2D and 3D simulations of semiconductor devices 
[18]. It has two components: ATHENA, which is a 2D 
process simulation framework with software tools for 
modeling semiconductor fabrication processes, and ATLAS, 
which predicts the electri- cal behavior of semiconductor 
structures and provides insight into the internal physical 
mechanisms associated with device operation. ATLAS 
produces outputs in three forms: the run time output as a 
guide to the progress of the simulations running, log files 
storing summaries of the electrical output information, and 
the solution files storing the data related to the values of 
the solution variables within the device. The numerous 
carrier transport models, the availability of differ- ent 
material substrates, and the ability to account for other 
physical effects (e.g. surface charges, lattice heating, 
genera- tion and recombination) have added greater 
flexibility to this device simulator. 

Semiconductor device operation is modeled in ATLAS by 
a set of one to six coupled, non-linear partial-differential 
equa- tions (PDE) [18]. ATLAS provides numerical solutions 
of these equations by calculating the values of the 
unknowns on a mesh of points within the device. The 
original, continuous model (similar to the PNP equations) is 
converted to a dis- crete, non-linear algebraic system by an 
internal discretization procedure. The sets of PDEs, the 
mesh and the discretiza- tion procedure determine the 
non-linear algebraic problem to be solved. This is solved 
using an iterative procedure that refines successive 
estimates of the solution. Different solution procedures 
(e.g. Gummel, Newton, or block iteration) exhibit different 
behavior with respect to convergence, accuracy, effi- 
ciency and robustness. 
2.3. KcsA ion channel simulations with SILVACO 
Based on its structure revealed by X-ray crystallography 
[2,13], we have built a KcsA channel model as shown in 
Fig. 1. The potassium channel is modeled here as a 
transmembrane lumen, surrounded by protein walls and 
with cylindrical reser- voirs of potassium and chloride ions 
at two ends [4-13]. The channel extends from 10 to 50 A ; 
with a narrow selectivity fil- ter of radius 1.5 A’ and length 
12 A’ and a wider segment of 
(1) and (2) iteratively and obtain self-consistent solutions for 
the potential, concentrations and fluxes of the ions along the 
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length 23 A’. The selectivity filter extends into the 
extracellular space, while the wider segment tapers 
towards the intracellu- lar space. Each reservoir at either 
end has a radius of 40 A’ and a width of around 8A’. 

To replicate the shape of the KcsA channel, the 
structure is made of two materials: one to mimic the 
water continuum that is conducting, and the other to 
mimic the non-conducting pro- tein walls [15,19]. The 
dielectric constants « of the conducting and _ non- 
conducting mediums are set as 80 and 2, respectively. 
The entire device is constructed from 55 rectangular 
regions, which is the maximum limit for ATLAS [18]. 
Each conduct- 
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Fig. 1 - The model of a KcsA ion channel built 
in SILVACO. The channel extends from 10 to 50 
A , with a narrow selectivity filter of radius 1.5 
A’ and length 12 A’ and a wider segment of 
length 23 A’. The selectivity filter extends into 
the extracellular space, while the wider 
segment tapers towards the intracellular 
space. 


ing region has the capability of being defined with its 
unique material properties. To this end, the carrier doping 
and their diffusion coefficients are specified in each 
conducting region as shown in Fig. 2. This allows us to 
emulate the fact that there are considerably fewer carriers 
in the vestibule of the channel compared to in the 
reservoirs. Interestingly, the selectivity fil- ter hardly has 
four to five ions at one instant [13], which is represented 
as an undoped region in Fig. 2. After the defini- tion of the 
device structure, the mesh points are specified to 
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Fig. 2 - The doping profile of the various 
conducting 

regions in the KcsA model. The extracellular 
baths have a symmetric concentration of 100 
mM equivalently, while the selectivity filter is 
left undoped. The profile to chosen to represent 


denote the locations where solutions are to be determined. 
The mesh is defined by a series of horizontal and vertical 
lines and the spacing between them. With a maximum 
limit of only 9600 mesh nodes in ATLAS [18], it is important 
to use the mesh lines wisely and in locations which are 
crucial for the device behavior. Two electrodes, the source 


and drain, each of length 20 A’ are placed at either ends 
of the reservoirs for applying a potential and measuring 
the electrodiffusion characteristics of the device. 

Despite the fact that ion channels can be modeled as 
nanoscale electronic devices as discussed above, there are 
sev- eral issues to be addressed in their electrical 
simulations. The carriers in electronic devices are electrons 
and holes, while those in ion channels are charged ions 
bound by hydration shells [10]. The masses and radii of the 
semiconductor carriers are significantly small compared to 
those of ions in channels. Furthermore, the diffusion 
coefficient of electrons is around 
500 cm? s“?, while that of K* ions is around 5 x 107’ cm? s7?. 
The energy band structure and crystalline structures of 
solid- state materials is well known, whereas that of liquids 
is less known. On the brighter side, the drift-diffusion PNP 
theory still holds for both solid-state devices and ion 
channels, providing a simplistic, macroscopic study of the 
flow of electrons and ions alike [15-17]. With this basis, we 
have incorporated many of the characteristics of ion flow 
into our channel model. Even though it is difficult take into 
account the large mass and radii of ions, we have adjusted 
the diffusion coefficients of elec- trons in each conducting 
region to describe the slower motion of ions in a channel’s 
vestibule [19]. Also, the net carrier doping and the degree 
of ionization is specified in a manner to pre- dict decreased 
concentrations in the channel’s vestibule and to obtain 
realistic current-voltage plots, as shown in Fig. 2. 

For this purpose, we have used the fact that the carrier 
doping density (per cm?) is given by n = 10? NayC, where 
Nay is the Avogadro’s number and Cis the concentration of 
the solution (mol/l). As will be demonstrated in the next 
section, with these modifications and appropriate mesh 
definition, it is possible 

to obtain realistic simulations results. 


3. Results 


The solid-state KcsA channel model was_ simulated 
using the PNP theory under various symmetric bath 
concentrations and_ potentials. The discretized PNP 
equations were solved by ATLAS using Gummel iteration 
and with each simulation run- time of <30 s. 

Figs. 3 and 4 show the potential contours along the 
chan- nel at a bath concentration of 100 mM and voltages 
of 100 and 200 mV, respectively. In Fig. 3, for example, a 
voltage potential of 100 mV was applied to the drain 
electrode in the reservoir at the right side, while the source 
electrode at the left side reservoir was kept at ground. The 
potential drops gradually along the channel, with most of 
the significant voltage drop 
occurring within the 12 length of narrow selectivity filter. 
A 


the realistic scenario in the vestibule of the 
channel [13]. 
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Figs. 5 and 6 plot the electric field contours along the 
channel axis at 100 and 200 mV, respectively, and a 100 
mM bath con- centration. The electric field variations can 
be observed along the channel, with noticeable changes 
at either mouths of the narrow selectivity filter. The 
entrance and exit of the channel 
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Fig. 3 - The SILVACO contour profile of the 
potential variation along the channel axis with 
100 mV applied across the channel and 
extracellular bath concentrations of 100 mM. 


are made rounded, as any sharp boundaries cause 
problems while solving the Poisson’s equation and lead to 
unwanted peaks. 

Figs. 7 and 8 plot of the electric fields and potentials, 
respec- tively, along the channel for varying drain voltages 
and a bath concentration of 100 mM. As observed, the 
simulation pro- gram takes into account an inherent built- 
in potential of the structure, which is around ~0.575 V in 


this case. The origin 
of this potential is similar to that of Nernst Potential in ion 


channels as the SILVACO program considers the same 
Boltz- mann energy distribution of carriers [10]. In this case, 
the bath concentration is 6 x 107° cm? (or 100 mM 
equivalently) and the inside of the channel has an 
intrinsic concentration of 
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Fig. 4 - The SILVACO contour profile of the 
potential variation along the channel axis with 
200 mV applied across the channel and 
extracellular bath concentrations of 100 mM. 
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Fig. 5 - The SILVACO contour profile of the 
electric field variation along the channel axis 
with 100 mV applied across the channel and 
extracellular bath concentrations of 100 mM. 


10*° cm-3. Plugging this in the formula for the built-in 
poten- tial V = (k7/q) In(Goath/Cchannei), WE get a value of 
~0.575 V. In Figs. 9 and 10, we plot the electric fields and 
potentials, respec- tively, along the channel for a fixed 
drain voltage of 100 mV and varying symmetric bath 
concentrations. As expected, the built-in potential 
increases with increasing extracellular con- centrations. 

In an attempt to check the validity of our model, we 
com- pare the current-voltage (IV) characteristics from our 
simula- tions with those of experiments [20]. Fig. 11 shows 
the corre- sponding IV plots of a KcsA channel at two 
symmetric bath concentrations of 250 and 500 mM. Using 
reverse engineer- ing, the material properties of the solid- 
state KcsA structure 
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Fig. 6 - The SILVACO contour profile of the 
electric field variation along the channel axis 
with 200 mV applied across the channel and 
extracellular bath concentrations of 100 mM. 
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Fig. 7 - Plots of electric field variation along 
the channel axis at varying voltages applied 
across the channel and extracellular bath 
concentrations of 100 mM. 


same 


have been adjusted so as to mimic the realistic ion-channel 
structure and obtain a good fit between the IV curves. As 
seen from Fig. 11, the simulated IV curves closely agree 
with the 


experimental IV curves until 70 mV. Above this voltage, 
there is a slight tendency towards saturation in the 
experimental curves, whereas the simulated IV curves 
keep rising almost linearly. Nevertheless, this degree of 
agreement between our simulated results and the 
experimental IV curves is much bet- ter than those reported 
in the previous literature [9]. 

Next, we consider the effect of surface charges on the 
electrical properties of KcsA channels. These surface 
charges are known to influence the gating, conductance, 
and toxin- binding effects of KcsA ion-channels, and their 
actual physical locations are presumed to be in the 
selectivity filter [21]. Being negatively charged, these 
surface charges aid the permeation of cations and impede 
the flow of anions through the KscA channel. Here, we 
explored the effects of such charges on ion permeation by 
placing three positive charges in the selectiv- ity filter at 
the interface of water and protein walls (at 40, 44, 
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Fig. 9 - Plots of electric field variation along the 
channel axis at varying extracellular bath 
concentrations and an applied voltage of 100 
mV across the channel. In each case, 


both the baths are symmetric, having the 


concentration. 


and 48 
A 
Fig. 8 - Plots of potential variation along the 
channel axis at varying voltages applied across 
the channel and extracellular bath 
concentrations of 100 mM. 


in Fig. 2). It is worth mentioning that the surface 
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charges in our simulations are positive (each Qc = 10” 


cm~*) 

since the majority carriers used in our model are 
electrons. The SILVACO program is versatile in allowing us 
to incorpo- rate any fixed charges, interface traps, or even 
bulk traps at the semiconductor-insulator interface and 
observing their effects on the conduction of a device [18]. 
Fig. 12 plots the poten- tial variation along the channel 
with a bath concentration of 100 mM and drain voltages 
of 0 V, 150 mV and 250 mV. In comparison to its steadily 
increasing nature in Fig. 7 (with no charges), the axial 
potential variation in Fig. 12 (with surface charges) has a 
plateau inside the selectivity filter, with three small 
humps resulting from the charges at those locations. The 
simulated IV curves (with and without charges) are plot- 
ted in Fig. 13 for a 250 mM bath concentration. Besides 
the 


Fig. 10 - Plots of potential variation along the 
channel axis at varying extracellular bath 
concentrations and an applied voltage of 100 
mV across the channel. In each case, both the 
baths are symmetric, having the same 
concentration. 
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Fig. 11 - The current-voltage characteristics of 
the KcsA channel model at two different 
symmetric bath concentrations (250 and 500 
mM). Comparison is made between the 
SILVACO simulated plots (—, 250 mM; - - -, 500 
mM) and the experimental data (*, KcsA: 250 
mM; O, KcsA: 500 mM). 


general increase in conductance, there is an added 
asymme- try in the linearity of the IV curves (in positive 
and negative voltage ranges) by surface charges. As seen 
in Fig. 13, the current component due to the majority 
carriers (electrons) is enhanced much more compared to 
that due to the minority carriers (holes) owing to the 
attractive and repulsive forces, respectively, exerted by 
the surface charges. 


4. Discussion 


We have demonstrated the possibility of representing a 
KcsA channel as a solid-state device, and employing a 
device simu- 
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Fig. 12 - Plots of potential variation along the 
channel axis at varying voltages applied 
across the channel, in the presence of surface 
charges, and extracellular bath concentrations 
of 100 mM. Three pairs of charges (each Qsc = 
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Fig. 13 - The simulated current-voltage 
characteristics of the KcsA channel model in 
the presence of surface charges and at two 
different symmetric bath concentrations (250 
and 500 mM). Besides the increase in 
conductance, there is an enhancement in the 
majority carrier current and decrease in the 
minority carrier current. The nature and 
amount of surface charges is same as that in 
Fig. 12. 


1014 cm~2) are located inside the selectivity 
filter at the interface between the water and 
protein walls (at 40, 44, and 48 A’ in Fig. 2). 
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lator, SILVACO, to simulate its electrical characteristics. 
With a goal to mimic the realistic and computationally 
feasible scenario inside an ion channel, we created a 
model with well-adjusted material parameters and 
employed SILVACO to obtain self-consistent solutions of 
the axial potential and ion fluxes. The simulated IV 
curves closely match the available experimental results, 
supporting the validity of the model. 

Besides giving insight into the general ion permeation 
mechanisms in a channel's vestibule, the simulation 
results do open exciting opportunities to extend this 
work. Addi- tional features can be incorporated in the 
model to bridge the gap between the nature of a 
semiconductor-insulator and a_ water-protein wall 
interface. In solid-state devices, the inter- face between 
a semiconductor and an overlying insulator is sharp and 
abrupt, while the dielectric boundary at a channel's 
water-protein wall interface is not so distinctly defined. 
This dielectric boundary and its curvature is important in 
decid- ing an ion’s potential energy, which varies as 
1/Ewater iN the 
baths and as 1/€protein inside the channel [22]. Even 
though the precise value Of €Eprotein is not crucial, it 
would better our 
model to put an additional protein region (eprotein ~= 10) 
at the 
dielectric boundary to smoothen the sharp water-protein 
wall interface. Also, it is suspected that the dielectric 
constant of water decreases inside the channel, and 
recently Monte Carlo simulations have shown that 
lowering its dielectric constant in a channel's vestibule 
alters a channel's selectivity for dif- ferent ions [5]. 
Studies can be performed to investigate the effects of 
changing €Ewate on the channel permeation and_ its 
selectivity for monovalent or divalent ions. 

Although our model is based on macroscopic drift- 
diffusion computations, they do_ describe the 
experimental IV data surprisingly well. It is, however, 
tempting to be able to perform simulations at molecular 
and atomistic levels (Monte Carlo and _ Brownian 
dynamics) with the computational ease of PNP programs 
[19]. In this context, a common platform can 
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be built which has the capability of running molecular-level 
simulations for a short timeframe, inferring the ion- 
transport parameters, and feeding them to the SILVACO 
program to cal- culate the ion fluxes within reasonable 
time. Such a hybrid model would blend the accurate 
predictions of molecular-level simulations with the 
computationally cheap macroscopic flux calculations. 
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